function r = RandNormBM( mu, sigma )
    while(1)
        v1 = 2*rand() - 1;
        v2 = 2*rand() - 1;
        s = v1^2 + v2^2;
        if(s<=1)
            break;
        end
    end
    
    r = sqrt((-2)*log(s)/s) * v1;
    
    r = r * sigma + mu;
end

